% Last modified by Heather Snell December 18, 2020.

%Take draws from the parameter space. 
%These draws are held in c which is size NVxNPxNDRAWS
%Then calculate the probability of each person's sequence of choices at
%each draw for that person. 
%These probabilities are held in PROBS which is size NPxNDRAWS

% Make sure rate is in NVth place
rrate = squeeze(c(:,NV,:));  %NP by NDRAWS
if discountfun=="QuasiHyper"
    QH = f;                  % Beta for hyperbolic
end

%computes the discount factor Psi

Psi = zeros(NDRAWS, NP, T+1);
if discountfun == "Exponential"
    for ii=0:T
        Psi(:,:,(ii+1))=((1./(1+rrate)).^(ii))';
     end
elseif discountfun == "HM"
    for ii=0:T
        Psi(:,:,(ii+1))= 1./(1+(rrate.*ii))';
    end
elseif discountfun == "Harvey"
    for ii= 0:T
        Psi(:,:,(ii+1))=(1/(1+ii)).^(rrate');
    end
%%% Quasi Hyper would take more set up work...
 elseif discountfun =="QuasiHyper"
    Psi(:,:,1) = 1;
    for ii=1:T
       Psi(:,:,(ii+1))= QH .* ((1 ./(1+rrate')).^(ii));
    end
else
   disp('Discount.fun not found') 
end
    

Psi = reshape(Psi, 1,1,NDRAWS,NP,T+1);

DX = zeros(NV-1,2,5,777,NDRAWS);

%%%% XX must have the data in it
for ii=1:NDRAWS
        DX(:,:,:,:,ii) = sum(bsxfun(@times, XX, Psi(1,1,ii,:,:)),5);
end



%Add up the non-price and non-rate variables times their WTP

%slightly different form for the ANA models
if (ANA==2 || ANA==3)
    c=permute(c,[2,1,3]);
    v=sum(bsxfun(@times, DX(1:(NV-1),:,:,:,:), reshape(c(1:(NV-1),:,:), NV-1,1,1,NP,NDRAWS)),1);      %sum(NROWSx(NV-1)xNDRAWS,2) gives NROWSx1xNDRAWS
else
    rscale = squeeze(c(:,1,:));
    c=permute(c,[2,1,3]);
    v=sum(bsxfun(@times, DX(2:8,:,:,:,:), reshape(c(2:8,:,:), 7,1,1,NP,NDRAWS)),1);      %sum(NROWSx(NV-1)xNDRAWS,2) gives NROWSx1xNDRAWS
end

v=squeeze(v);                              %NROWSxNDRAWS
if (ANA~=2 && ANA~=3)
    v=v-squeeze(DX(1,:,:,:,:)); %Subtract price
    v=bsxfun(@times, reshape(rscale, 1, 1, NP, NDRAWS),v) ; %Multiply everything by price/scale coef
end

v=exp(v); 
v(isinf(v))=10.^18; 
clear Psi rrate DX


